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ABSTRACT 

We  consider  wave  equations  with  damping  in  the  boundary  conditions. 
Techniques  to  ascertain  the  uniform  preservation  under  approximation  of  ex¬ 
ponential  stability  are  presented.  Several  schemes  for  which  preservation  can 
be  guaranteed  are  analyzed.  Numerical  results  that  demonstrate  the  lack  of 
stability  under  approximation  for  several  popular  schemes  (including  standard 
finite  difference  and  finite  element  schemes)  are  given. 
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1  Introduction 


In  this  paper,  we  consider  approximation  methods  for  the  boundary  damped  nor¬ 
malized  (i.e. ,  wave  speed  c  =  l)  wave  equation  system 


d?u(t,x)  =  A  zu(t,x), 

t  >  0, 

*  £  n  c  Rn, 

u(t,x )  =  0, 

t  >  0, 

x  £  r0  c  so, 

dnu(t,x)  +  adtu(t,x)  =  0, 

t  >  0, 

x  £  Tj  C  30, 

where  the  domain  O  is  an  open  bounded  subset  of  Rn,  To  is  a  relatively  closed  subset 
of  the  boundary  30  of  the  domain  0,  and  IT  is  the  complementary  subset  of  T0  in 
the  boundary  30.  The  symbol  3„  represents  the  directional  derivative  operator  in 
the  outward  normal  direction  of  the  boundary.  Since  the  boundary  condition  (1-2) 
is  often  referred  to  as  a  reflecting  boundary  condition,  To  is  called  the  reflecting 
boundary.  Similarly,  IT  is  called  a  partially  absorbing  boundary. 

The  system  (l.l)-(1.3)  arises  in  many  important  models  for  distributed  pa¬ 
rameter  control  problems.  In  particular,  in  the  model  of  a  vibrating  flexible  mem¬ 
brane,  the  solution  u(t,  x)  represents  the  transverse  displacement  of  the  membrane, 
and  in  models  for  acoustic  pressure  fields,  the  solution  u(t,z)  represents  the  fluid 
pressure  (see,  [BKS],  [BKSW],  (BPS],  [Li],  for  more  detailed  examples).  It  can 
be  shown  that  for  a  given  initial  state  u(0,x)  =  u o(z)  and  ut(0,x)  —  u0(z),  the 
solution  of  (1.1)— (1.3)  decays  exponentially  in  time  and  the  decay  rate  is  uniform 
for  all  initial  states  (u0,v0)  in  a  certain  function  space  (see  [C],  [L]).  The  stability 
of  the  solutions  of  (l.l)-(l.3)  plays  an  essential  role  in  several  control  theoretical 
issues  ([BK],[BW]).  In  this  paper,  we  are  interested  in  approximation  methods  that 
uniformly  preserve  the  exponential  stability  of  the  solutions  of  (l.l)-(1.3)  for  the 
approximate  solutions. 

The  equations  (l.l)-(l.3)  can  be  written  abstractly  as  a  differential  equation 

~w(t)  =  Aw(t),  t  >  0,ie(t)  £  )( , 

at 

in  an  infinite  dimensional  function  space  K.  In  this  context,  we  consider  a  linear 
control  system  given  by 

(1.4)  ~w{t)  =  Aw{t)+Bh{t),  h(t)  £  Rm, 

where  h  is  a  control  input  and  B  is  a  linear  operator  from  Rm  into  "H .  The  most 
common  approach  for  the  approximation  of  a  control  problem  involving  (1-4)  is  to 
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formulate  a  sequence  of  finite  dimensional  control  systems  of  the  form 

(1.5)  wN{t)  =  A NwN(t)  +  BNk{t),  t>  0,  wN{t)eUN, 

where  the  dimension  of  the  space  )lM  increases  toward  infinity  as  N  tends  to  in¬ 
finity.  In  general,  equation  (1.5)  is  derived  from  (1.4)  using  space  discretization 
techniques  such  as  finite  difference,  finite  elements  or  spectral  methods  developed 
for  the  approximation  of  the  solutions  of  (l.l)— (1.3).  A  control  strategy  is  then 
designed  for  the  finite  dimensional  control  problem  involving  (1.5).  This  control 
is  used  as  an  approximation  to  the  desired  control  function  for  the  control  prob¬ 
lem  involving  (1.4)  (for  example,  see  [Gj,  [BK],  [BWj).  One  of  the  most  practical 
conditions  to  assure  the  weil-posedness  of  the  finite  dimensional  control  problems, 
as  well  as  the  convergence  of  the  approximate  controls  to  the  desired  control  for 
the  infinite-dimensional  system  is  that  the  solutions  of  (1.5)  for  k  =  0  preserve  the 
exponential  decay  of  the  solutions  of  (1.1)-(1.3). 

From  a  stability  analysis  of  the  solutions  of  (l.l)  —  (1.3),  it  is  easy  to  see  that 
the  energy  dissipation  in  (l.l)— (1.3)  comes  exclusively  from  the  boundary  condition 
(1.3).  Since  no  medium  damping  exists  in  (l.l),  we  refer  to  (l.l)— (1.3)  as  a  weakly 
damped  wave  equation.  Although  many  approximation  techniques  can  provide  con¬ 
vergent  approximations  for  the  solutions  of  (l.l)— (1.3)  by  the  solutions  of  (1.5)  with 
h  =  0,  the  nature  of  the  dissipation  in  (l.l)— (l .3)  makes  it  very  difficult  to  preserve 
the  uniform  exponential  decay  of  the  solutions  of  (l.l)  — (1.3).  Numerical  results 
from  our  investigations  reveal  that  most  of  the  popular  discretization  techniques 
can  not  maintain  a  uniform  decay  rate  in  the  solutions  of  (1.5)  with  h  ~  0  as  the 
dimension  of  the  approximate  system  (1.5)  increases. 

Although  the  preservation  of  the  stability  or  the  stabilizability  of  hyperbolic 
type  control  systems  is  wrel!-known  to  be  a  delicate  approximation  problem,  there 
exist,  to  our  knowledge,  only  a  few  analytical  results  and  these  are  for  the  approx¬ 
imation  of  “hyperbolic”  delay  equations  (see  [IK]).  In  this  paper,  we  present  a 
general  approach  for  the  analysis  of  uniform  preservation  of  exponential  stability  of 
approximation  systems  for  weakly  damped  wave  equations.  This  approach  is  later 
used  to  show  that  a  particular  mixed  finite  element  method  and  ;  dynomial  based 
Galerkin  methods  preserve  a  uniform  exponential  decay  rate  in  the  approximate 
solutions  as  the  dimension  of  the  approximate  system  increases  to  infinity. 

An  outline  of  the  paper  is  as  follows.  In  Section  2,  we  specify  the  class  of 
approximation  methods  for  the  weakly  damped  wave  equation  considered  in  this 
paper.  Then  wc  present  a  general  approach  for  the  analysis  of  decay  rate  in  the 
solutions  of  these  approximate  wave  equations.  In  Section  3,  the  approach  developed 
in  the  proceeding  section  is  used  to  prove  that  a  mixed  finite  element  method  for 


the  1-dimensional  problem  preserves  uniformly  an  exponential  decay  rate  in  the 
approximate  solutions.  In  Section  4,  polynomial  based  Galerkin  approximations  of 
weakly  damped  wave  equations  in  hypercubical  domains  are  analyzed.  These  are 
also  shown  to  preserve  uniformly  an  exponential  decay  rate.  In  Section  5,  several 
popular  discretization  methods  are  investigated  numerically;  a  sharp  distinction 
between  the  methods  analyzed  in  the  earlier  sections  and  the  other  popular  methods 
emerges  from  our  numerical  findings.  In  Section  6,  our  concluding  remarks  offer  a 
perspective  of  the  results  presented  here. 

We  finish  this  introduction  by  observing  that  interest  in  the  weakly  damped 
wave  equation  is  also  strongly  motivated  by  the  question  of  exact  controllability  of 
the  wave  equation  via  partial  boundary  control.  As  it  was  pointed  out  by  Datko  (see 
D),  exact  controllability  or  stabilizability  via  boundary  control  may  be  extremely 
sensitive  to  perturbations.  However,  we  view  equations  (l.l)-(l .3)  as  the  model 
of  a  physical  system,  where  the  control  input  is  a  nonhomogeneous  term  in  (1.1) 
as  it  is  formulated  in  ([BFSj,  [BKSW],  [BKS],  etc.).  Thus,  our  control  input  does 
not  introduce  changes  in  the  boundary  condition  (1.3).  More  discussion  on  the 
significance  of  our  results  for  the  preservation  of  boundary  stabilizability  of  the 
wave  equation  is  given  in  Section  6. 
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2  Estimation  of  decay  rate  for  approximate  solutions  of  the 
wave  equation:  A  general  approach 

In  many  analytical  studies  of  control  problems,  the  system  (l.l)-(l .3)  is  taken  in 
the  sense  of  distributions  and  one  seeks  mild  or  weak  solutions.  As  a  second  order 
equation,  the  state  space  for  the  mild  solution  (u(<),u,(0)  °f  (l-l)-(l-3)  is  taken  to 
be  )l  —  Hpa{U)  x  L2(n),  where  the  Hilbert  space  H^g(Cl)  is  defined  by 

Hr0(n)  =  MO  €  ff1(n)|ti(i)  =  0,x  <E  r0}. 

For  simplicity  of  analysis,  we  assume  that  r0  has  a  non  empty  relative  interior;  thus 
an  inner  product  in  H^0(Cl)  that  is  equivalent  to  the  usual  H1( fl)  inner  product  can 
be  defined  by 

<u,v>!=  J  Vu(i)  •  Vv[x)dx,  u,v  €  jFfro(fl). 

The  infinitesimal  generator  for  the  mild  solution  semigroup  of  (l.l)-(l.3)  is  defined 
by 

D!A)  =  I  (“M€k|u€i^0(fl),Au6L2(n),  1 
{  dnu{x )  -f  av(x )  =  0,  x  6  T2  J 

A  :  D(A)  C  M  •  *■  1/ ,  A(u,v)  =  (u,  Axu). 

It  is  easy  to  verify  that  A  generates  a  Co-contraction  semigroup  5(i)  in  'A.  Ii  is 
shown  in  [C]  and  [L]  that  there  exist  constants  M  >  l,w  >  0  such  that 

||5(t)|U(»)  <  Me'wl,  for  t  >  0. 

Alternatively  and  equivalently,  for  the  initial  state  (u0,v0)  €  V{A),  the  solution 
(u(t),v(t))  =  5(f)(u0,u0)  satisfies  the  following  variational  equality, 

(2.1) 

for  all  ( f,g )  €  H^o(Q)  x  H^n[ fl).  Equation  (2.1)  is  often  referred  to  as  a  weak 
formulation  of  ( 1  - 1)-(  1 .3)  and  it  is  a  natural  formulation  in  which  to  consider  ap¬ 
proximation  methods. 

Consider  two  finite  dimensional  spaces  and  H j  of  functions  defined  on 
0  given  by 

=  sPan{^}JL|f  H?  =  spanfO^j. 


<«(*)./>  1  =  <v(t),f> ! 

(*h?>£J(n)  =  ~<u{t),g> \-j  av(f,x)p(x)a!x 

i 
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We  consider  approximations  to  the  solutions  of  (1.1)-(1.3)  in  the  form 


N 


N 


N 


*  =  Y2uk[*)<f>k{x)>  vNi^x)  =  YlvkW'Pki*)- 

Jt  =  l  k=l 


In  most  cases,  and  are  chosen  to  be  subspaces  of  //po(Q)  and  L2(0), 
respectively.  However,  this  condition  is  not  essential  in  the  analysis  presented  below. 
Let  us  denote  the  column  vectors  of  the  coefficients  (<)  by 


( <(<) ) 

(  v?(t)  \ 

Si 

* 

II 

II 

■vu 

* 

1 

V  > 

t  *#(0  ) 

Then  system  (l.l)-(l.3)  is  approximated  by  an  ordinary  differential  equation  of  the 
form 


(2.2) 


'  KN 

0 

(  *"(0  ]  _  ' 

0 

MN  _ 

l  j 

0 

-kn 


Kn 

~Bn 


*N(t)  ) 

vN(t)  )  ' 


The  matices  KN  and  MN  are  symmetric  and  positive  definite  by  construction. 
The  matrix  BN  is  symmetric  nonnegative  definite.  It  arises  from  the  boundary 
integral  term  in  (2.1).  In  general,  mappings  qN  :  H *  >—>  are  chosen  [IK].  The 
approximate  solution  (us(t),  vN(t))  is  required  to  satisfy  the  following  variational 
system  which  is  analoguous  to  (2.1): 


!<> .»((),<  >IW  =  /r^"(vv(())9',«)ix. 

In  particular,  if  the  mappings  qN  are  defined  as  qN(ip^)  —  cf>^ ,  j  =  1 N,  then, 
the  matrices  KN ,M'V  and  BN  are  defined  by 


K”  =  /nV^f(x)-V^(x)dx, 

M"  =  Jn  [x)dx, 

B"  =  Jr  a<f>? {x)(t>f  (x)dx. 

In  the  case  of  Galerkin  type  of  methods,  i.e.,  finite  elements,  spline  based  Galerkin 
methods,  polynomial  based  Galerkin  methods,  =  //*,  and  therefore,  qK  —  I 
the  identity  operator. 
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We  are  interested  in  the  exponential  decay  rate  of  the  approximate  solutions 
( uN ,vN ).  First,  we  need  to  assume  that  the  appropriate  norm  for  the  approxi¬ 
mate  solutions  is  used.  In  general,  the  convergence  properties  of  the  approximation 
schemes  can  be  stated  as  follows:  Let  iN  be  a  mapping  from  =  H*  x  H ?  into 
U.  Then  (uN,vs)  is  said  to  converge  to  (u,v),  if 

as  N  tends  to  infinity. 

Definition  2.1  (Uniform  exponentially  stable  approximation )  A  given  approxima¬ 
tion  method  is  said  to  preserve  uniformly  the  exponential  stability  of  the  solutions 
of  (l.l)-(l.3),  if  there  exist  constants  M  and  a  >  0  independent  of  N  such  that  for 
any  initial  state  [u0,va)  the  corresponding  aproximate  solutions  satisfy 

\'iN{uN(t),vN(t))h  <Me-af||(u0,v0)i|.v,,V  =  1,2,.... 


It  follows  from  the  positivity  of  the  matrices  KN  and  MN  that  a  norm  is 
defined  on  Hf  x  H F  by 

,!(uA  ,VN)\\hn  xHn  =<  KNUN,UN  >  rn  +  <  MNVS ,  V'  >ftN, 

where  us  and  are  vector  representations  of  uN  and  vN  with  respect  to  the  bases 
of  //{v  and  H? ,  respectively.  We  require  the  following  condition  to  hold. 

Condition  2.1  (Norm  compatibility)  There  exist  constants  ci,c2  independent  of  N 
such  that  for  all  {uN ,vN)  6  x  H? , 


j(u 


.V  Va 


iHf  xH"  — 


<  II*  (w  ,  v  )  [|x  <  c2||(«  ,  ) 


i  ■HfxH?- 


We  note  that  in  the  case  of  Galerkin  type  of  methods  mentioned  previously, 
the  above  condition  holds  with  cx  =  c2  =  1.  Under  the  above  condition,  in  order 
to  establish  preservation  of  exponential  stability,  it  is  sufficient  to  show  that  the 
approximate  solutions  (us(t),  vN  (t))  have  a  uniform  exponential  decay  rate  under 
the  norm  ]j  • 

In  the  analysis  of  the  stability  of  solutions  of  ( 1  ,l)-(  1 .3) ,  a  generalized  Lyaponov 
function  Q(t)  is  used  in  ([Cl,  [ L]) .  It  is  defined  by 

0(0  =  ^(’:(ti(0,u(0)!iv)  +  /n  |2v(t) (£  •  Vtt(i))  +  u(Ou(o|  dx 
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where  £  is  a  vector  field  defined  on  ft  that  represents  certain  characteristics  of  the 
domain  ft  and  its  boundary,  £,-_,- (x)  =  dZi£,(x),  and  the  term  E(t)  =  ||(u(t), u(t)j| \ 
represents  the  energy  associated  with  a  solution  ( u(t),v(t ))  at  time  t.  The  main 
steps  of  the  analysis  consist  in  establishing  the  following  results: 

(LI)  There  exists  a  constant  Tx  >  0  such  that  for  all  t  >  Tx, 

o(o  >  \m- 

(L2)  There  exists  a  constant  T2,  such  that  for  all  t  >  T2,  Q(t)  <  0. 

(L3)  For  any  given  T  >  0,  there  exists  a  constant  M  such  that  for  all  t  <  T, 
\Q(t)\<ME(t). 

From  (Ll)-(L3),  we  can  find  a  constant  C  such  that  for  all  initial  conditions  and 
for  all  t  >  T  =  m&x{Tx,T2}  we  have 

£(1)<1W)<cM 

W  ~  t  ~  t 

Using  the  semigroup  property  of  5(t),  e.g.,  see  [P,  p.116],  we  conclude  that  S(t )  is 
exponentially  stable  as  t  tends  to  infinity. 

Since  (2.2)  is  an  approximation  of  (l.l)-(l.3),  it  is  natural  to  attempt  to 
follow  a  similar  idea  to  analyze  the  stability  of  the  solutions  of  (2.2).  Thus,  we 
introduce  a  function  QN{t)  defined  by 

QN(t)  =  ±(<KV(t),u*(t)>RN  +<MV(1),^(1)  >rN) 

4-  <  WNuN(t),u*(0  >RN 

where  is  an  iV  x  jV  matrix.  Obviously,  WN  plays  a  role  analoguous  to  that  of 
the  integral  term  in  Q(t)  involving  the  vector  field  £(•).  Then  if  we  take  EN(t)  = 
j|(uw(t),  u/v(<))||J/N>;//N,  by  establishing  (Ll)-(L3)  for  some  constants  TX,T2,  and  M, 
independent  of  N ,  we  can  obtain  the  uniform  exponential  stability  of  the  approx¬ 
imate  solutions  following  an  argument  similar  to  that  outlined  previously  for  the 
solutions  of  (l.l)-(l .3).  The  following  condition  is  required  if  (Ll)  and  (L3)  hold 
with  constants  Tx  and  M  independent  of  .V. 

Condition  2.2  (Boundedness  ofWN)  There  exist  constants  (3X  and  Bi  independent 
of  N  such  that  that  for  all  fN ,glW  £  RN , 

|  <  WNfN,gN  >RN  |</?,<  KNgN,gN  >rN  +&  <  MN  fN ,  fN  >R„  . 


7 


Lemma  2.1  Under  Condition  2.2,  there  exists  a  constant  7\  independent  oj  N  such 
that  for  all  t  >T\, 

Q"(t)  >  1). 

Moreover,  for  any  given  T  >  0,  there  exists  a  constant  M  independent  of  N  such 
that  for  all  t  <  T , 

|cj"(l)l  £  MEN(t). 

Proof:  It  is  sufficient  to  take  T)  =  4max{/?i,/?2}  and  M  >T/2  +  ma x{/?i,/52}- 

□ 

As  it  is  in  the  case  of  the  analysis  of  the  decay  rate  for  solutions  of  (l.l)-(1.3), 
condition  (L2)  is  the  most  difficult  to  establish  (see  [C],  [L]).  In  fact,  by  using  (2.2), 
we  can  compute  the  derivative  of  QN  as  follows 

$*(«)  =  A(<  Knun,  un  >rn  +  <  MNvN ,  vN  >rn) 

-t  <  Bnvn,vn  >rs  +  <  WNvN,vN  >RN 
-  <  WS{MN)~\KNuN  -r  BNvN),uN  >r»  . 

The  above  derivation  involves  use  of  the  expressions 

hN(t)  =  vN(t),  iN{t)  =  -{MN)~1(KNuN(t)  +  BNvN(t)). 

obtained  from  (2.2).  By  adding  and  subtracting 

(<  K” uN,c"  >„ »  +  <  Mnvn,vn  >rn)J2, 

we  obtain 

(2.3)  QN{t )  =  KNuN{t),us{t)>Rs  +  <  MNvN{t),vs{t)  >rN) 

+  <  KNuN[t),uN{t)  >rN  -  <  KNuN[t),  CNuN{t)  >Rs 
+  <  M"vs{t),vs{t)  >R»  +  <  MNvN{t),CNvN{t)  >Rs 
-t  <  BNvN{t),vN{t)  >R»  -  <  BNvN{t),CNuN{t)  >R», 

where  the  matrix  CN  is  defined  by  ( CN)T  =  Ws 

The  following  conditions  can  be  of  practical  use  in  establishing  (L2). 

Condition  2.3  There  exists  constants  a\,ai,az,p,  and  6,  with  (p  +  S)  <  j,  each 
independent  of  N  such  that 
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(SI)  For  any  fN  G  RN 


>„»  <  (a,  + 

+P  <  Mn fN ,fN  >rn; 


(S2)  For  any  fN,gN  G  i2  , 


<  JfY,,#  >*»  -  <  KNgN,CNgN  >R„  -  <  B»fN,CNgN  >RN 

<  (S  +  m)  <  >*.v  +-  <  BNfN,fN  >RN  . 

M 


We  note  that  a  necessary  (but  not  suffir  r„)  condition  that  (S2)  hold  (take 
(/N  =  0)  is 

(S3)  For  any  gN  G  RN , 

<  K  g  ,g  >rn  —  <  K  g  ,G  <7  >R"<  o  <  K  g  ,g  >R»  . 

If  a  matrix  VF^  can  be  found  such  that  Condition  2.3  holds,  we  can  show 
that  (L2)  holds  immediately. 

Lemma  2.2  Under  Condition  2.8  there  exists  a  constant  T2  independent  of  N  such 
that  QN(t)  <  0  for  all  t  >  T2. 

Proof:  Let  p  be  chosen  such  that  6  +  p  <  1/2,  then  by  taking  T2  =  a x  +  a2/ p  +  a3/p, 
using  (2.3),  we  have  QN  (t)  <  0  for  all  £  >  T2. 

□ 

In  order  to  obtain  a  uniform  exponential  decay  rate  for  the  approximate 
solution  semigroups  S^(-),  where  (uN(t),  v^t))  =  5Ar(t)(uA  (0),  vA  (0)),  we  also 
need  the  following  Lax-stability  condition  for  the  approximation  scheme. 

Condition  2.4  There  exist  constants  C  >  1  and  u>  >  0  independent  of  N  such  that 

<  CeWt 

for  all  t  >  0. 


The  above  condition  is  automatically  satisfied  if  the  approximation  is  con¬ 
vergent.  Finally,  we  state  the  main  result  of  this  section. 


9 


Theorem  2.1  Suppose  that  a  matrix  can  be  found  for  each  N  such  that  Con¬ 
ditions  2. 2-2. 4  hold.  Then,  there  exists  constants  L  >  1  and  a  >  0  such  that 


ll*5'N(0llt(K1NxH2A')  ^  Le  at 

for  all  t  >  0  and  for  all  N . 

Proof:  By  Lemma  2.1  and  2.2,  we  can  find  constants  T  and  M  independent  of  N 
such  that  for  all  t  >  T,  the  following  statements  hold. 

(a)  QN{t)>tEN(t)j4. 

(b)  QN(t)  <0. 

(c)  Q"{T)  <  MEN{T). 

By  Condition  2.4  and  (c),  we  obtain 

Qn{T)  <  MES{T )  <  CtMe2uTEy{ 0). 

Combining  the  above  inequality  with  (a)  and  (b),  we  obtain 

EN(t)  <  jC\Me2uTEN{ 0). 

If  we  let  A2  =  4 C‘Me2ujT ,  the  above  inequality  can  be  rewritten  as 


r(0(«*(o), *>))&»**»  <  Tl!(uiV(o),vAr(o))!|Jrf,xWN, 


for  all  (u'v(0), vw(0))  £  x  Hj  and  t  >  T.  As  a  consequence,  for  r  =  4A2,  we 
have 

for  all  N.  Using  the  semigroup  property  of  S^(-)  and  Condition  2.4,  we  obtain  the 
inequality 

11^(0 IL(h,nxh^’)  ^  L  '  e 

by  taking  a  =  In  2/4 A2  and  L  is  a  constant  independent  of  N. 

□ 

It  is  clear  now  that  the  difficult  part  of  the  analysis  for  a  given  approximation 
is  to  estabish  Conditions  2. 1-2.4.  In  the  sub  equence  sections,  we  will  specify  the 
required  matrix  IV'v  for  two  different  types  of  approximation  methods  and  verify 
i-tiat  Conditions  2. 1-2.4  hold. 


3  Mixed  finite  element  methods 


In  the  most  common  implementations  of  finite  element  methods,  the  two  compo¬ 
nents  of  the  mild  solutions  of  (l.l)-(1.3)  are  approximated  by  functions  uN((,x) 
and  vN(t,x)  with  the  same  smoothness  in  spatial  variable  x,  or  more  precisely,  the 
approximation  spaces  H jA  and  are  often  chosen  to  be  identical.  However,  both 

the  analysis  for  the  existence  of  the  solutions  of  ( 1 . 1  )-(l .3)  and  nature  of  the  wave 
propagation  suggest  that  u  and  v  have  different  smoothness  in  x.  The  methods 
analyzed  here  use  different  approximation  spaces  H ±  and  H?  which  give  differ¬ 
ent  smoothness  in  i  to  the  functions  uN  and  vN .  The  term  “mixed  finite  element 
method”  is  used  heie  to  indicate  that  two  different  types  of  approximation  elements 
are  involved  in  these  schemes. 

We  first  consider  the  case  of  a  one-dimensional  wave  equation.  Let  fi  = 
(0, 1),  so  that  the  equations  (l.l)-(1.3)  can  be  written  for  a  typical  set  of  boundary 
conditions  (i.e.  L0  —  {0},Fi  -  { 1 } ) 

(3-1)  —u(t,x)  =  t  >  0, x  6  (0,1), 

(3.2)  u{t,0)  =  0,  t>  0, 

d  d 

(3.3)  ^u(t,l)  =  -a— «(M),  t  >  0. 


The  following  approximation  method  has  been  proposed  by  Ito  and  Kappel 
in  [IK 1 .  The  basis  elements  6%  and  for  spaces  H?  and  77^,  respectively  are 
defined  by 


1  —  N\x  —  x*|,  x  6  [x*-i,x*+I]nS0, 1], 
0,  otherwise, 


and 


f  1,  x  6  (it_1,xA  +  1]  fli!0,  I], 
[  0,  otherwise  , 


where  k  =  1 ,  -  *  * ,  TV  and  x*  —  k/N.  It  is  easy  to  see  that  forms  a  basis 

for  the  set  of  continuous  piecewise  linear  polynomials  (i.e.  linear  splines)  on  [0,1] 
corresponding  to  the  uniform  mesh  {k/N}k= 0  with  function  value  equal  to  zero 
at  x  =  0.  The  s  are  piecewise  constant  functions  with  the  same  support  as 
the  corresponding  <+>%' s-  Therefore,  7/ C  7/£(f2)  =  {<p  €  7/ 1  (0, 1 ) |0(O)  =  0}  and 
77^  =  span{t/>*r}f=,  C  L2(fi).  The  matrices  MS,KN,  and  BN  of  the  previous 
section  are  defined  by 


[  ilf  {x)ip* (x)dx, 
Jo 


II 


ities. 


* s  =  £  i#v£*ii*)** 

B"  =  a*f(l)*"(l). 

Before  choosing  the  matrix  WN ,  we  first  establish  the  following  useful  equal- 


Lemma  3.1  For  any  fN  6  RN ,  the  following  equalities  hold. 

(i)  Let  gN  —  MN  fN ,  then  the  component  is  given  by 

I  (2/."  +  /,")/K  k  =  1 

9?  =  (/£.,  +  2/,"  +  /£,)/*,  k  =  2,...,N~l, 

U/»+/"-,)W  *  =  JV. 


fiii 


Let  —  KN fN ,  then  the  component  g £  is  given  by 
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\  N{2 /?-/?),  k  ~  1 

■  #(-/£, +2/4" -/&.),  k  =  2,...,N-l, 
NUR-fR..),  k  =  N. 


(iii)  For  any  fN  6  RN ,  we  have 


<  >„„  =  If/,")2  +  1  £(/f  + 

1  *= 2 


AT’ 


(iv)  For  any  fN  6  toe  have 

<  >*»  =  ^(/r)2 + n  E(/kN  -  fw. 

k  =  2 


Proof:  Equalities  (i)  and  (ii)  follow  directly  from  the  definitions  of  the  matrices  MN 
and  Kn ,  respectively.  The  following  equalities  are  used  frequently  in  this  section, 
and  in  particular,  are  useful  in  establishing  (iii)  and  (iv) .  Let  a*,  bk  be  given  numbers 


for  k  =  1, . . 

..,N.  Then 

TV- 1 

/V 

(3.4) 

a2bx  +  a^bN  +  X  (a*  +  ak+\)bk 

—  X  Ok(bk  +  bk~i), 

k  =  2 

k=  2 

n-  1 

TV 

(3.5) 

—  a26j  +  ajvbN  +  X(a*  —  ak+i)^k 

=  X  a*(&*  ~ 

k  =  2 

k  =  2 
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From  (i),  we  have 

<mn/',,/">s»  =  ^(2/,"  +  /")/"  +  ~ua  +  /;_,)/; 

+~  zW  +  +  (/£.  +  C))C- 

JV  k=2 

By  taking  ak  —  fj?  -f  fk_1,bk  =  ff,  and  applying  (3.4),  we  obtain 

<  M"f,f  >*„=  i(/f )’  +  i  E(/f  + 

iV  *  =  2 

From  (ii),  we  have 

<  KNfN,  fN  >*„  =  JV(2/,"  -  /«)/,«  +  N(fff  -  /£_,)/£ 

+N  £((/."  -  /?_,)  -  (/"  ,  - 

*=  2 

By  taking  ak  =  —  }k_x,  and  =  /^,  and  applying  (3.5),  we  obtain 

<  KNfH,fK  >*»=  N(fX)2  +  N  •£(/?  -  /£,)>. 

k-2 


□ 


Now  we  can  define  a  matrix  (and  hence  the  matrix  WN  =  ( CN)TM A) 
as  follows.  For  any  vector  fN  €  /2‘v,  the  vector  =  CN f  N  is  given  by 

[  fi,  k  =  1, 

<?"  =  *(/&»“ /£i).  *  =  2,...,JV-1, 

l  2Af(/£  -  /#_,),  fc  =  N. 


The  definition  of  can  be  seen  as  a  discrete  approximation  of  the  term  in  the 
multiplier  Q(t)  involving  the  vector  field  t.  In  fact,  in  our  case,  we  ca;.  take 
£(x)  =  2z.  In  fact,  for  any  function  fN{x)  =  T,k=i  fk  (x)  6  Hk,  we  approxi¬ 
mate  2 t(x)dfN{x)/dx  by 

eN(x)  =xPN{x)T -CNfN. 

where  r/jN (x)T  is  a  row  vector  valued  function  given  by  tpN{x)T  =  (xp^(x),. . . ,  rp%(x)). 
One  can  verify  from  the  definition  of  CN ,  4>k  and  <f>k  ,  that  on  each  sub-interval 
[x*-i>  x*],  for  k  =  3, . . . ,  N  —  1,  we  have 

fN{Xk)  -  fN{xk-i)  ,  /^(xjt+i)  -  /N(x*_i) 

Z*-l - (-  X* - 

%k  %k-2  %k+ 1  1 


nN 


(*) 
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which  can  be  regarded  as  a  discrete  approximation  to  the  first  integral  in  the  defi¬ 
nition  of  <?(£)• 

Thus,  for  any  vN  €  H?  with 

vN{x)  =  J2vk'l’kix)^ 

fc-i 

and  eN  as  defined  above,  the  following  equality  holds 

f1  eN(x)vN{x)dx  =<  WNvN,fN  >RN  . 

Jo 

The  Lax-stability  for  this  approximation  scheme  is  given  by  the  convergence 
analysis  in  [IK],  and  the  norm  compatibility  condition  holds  with  Ci  =  c2  ~  1.  We 
need  to  establish  Conditions  2.2  and  2.3. 

Lemma  3.2  For  any  vector  fN,gN  €  RN ,  the  following  inequality  holds 

\<MNfN,CNgN  >rn  |  <  2  <  MNfN,fN  >RN  +2  <  KNgN,gN  >RN  . 

Proof:  By  definition  of  the  matrices  and  CN ,  we  have 

<  M'7",cv  >„„=  ^(2/*  +  ff)if  +  2UH  +  fZ-M  -  9?-,) 

+4  Ek{(fM  +  / k)  +  (A*  -  ((9t+i  -  Ok)  +  [Of  ~  of-,))  ■ 

iV  k= 2 

Using  the  inequality  ja6|  <  ( a2/N  +  A*62)/2,  we  obtain 

(3.6)  [(gi  +  02) (b\  +  62)|  <  — (a2  +  a\)  +  N (62  T  h\). 

By  (3.6),  we  obtain 

|  <  MNfN,CNgN  >„»  |  <  ^  ((/IT  H/.'Wj'T  +  U" +  /»-,)’) 

£  (</&,  +  AT  +  (/,"  +  /,' 

+w  ((sD2  +  («r  -  jJT  +  (si?  -  9?-,)1) 

+JV  £  ((s&.  -  O’  +  (dk  -  PC'-,)2) 

Jt=2 
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<  |  (<a")2  +  E<a" + am2) 

+2A^  f(ff?)2  +  “  ffit^x)2) 

\  *=2 

Lemma  3.3  For  any  fN  £  the  following  equalities  hold. 

(3.7)  <  MNfN,CNfN  >RN  =  -<MNfN,fN>RN 

+^fs+fN-i)fs-{fN+fN-,)\ 

(3.8)  <KNfN,CNfN>R»  =  <KNfN,fN  >R, 


Proof:  First  consider  (3.7)  and  observe,  as  in  the  previous  proof,  that  we  have 

<  M* fN ,CN fN  >*„=  -t(2/«  +  +  2(/y  +  /#_,)(/#  - 

4  E  *  ((At,  +  A")  +  (A"  +  AM)  ((/»+,  +  A")  -  (A"  +  A"-,))  • 
JV  *  =  2 

The  sum  SXV  on  the  right  hand  side  can  be  written  as 

i  E  [(*  +  D(At,  +  A")2  -  MA"  +  AM2  -  (/it i  +  A")2] 

JV  i=2 

=  (fX  +  AM2  -  |<A"  +  A")2  -  £  (/*  +  /",)=. 

Therefore,  we  obtain 


<  m''/",cMv  >„„=  -—{fXV-jj  Et/t+AM'+'iUK +/*-,)/«-(/£+ AM2- 

Thus  using  Lemma  3.1  (iii),  we  find 

<  mn  fN ,  cNr  >RN=-<  MNr,  r  >RN  +4(/^r + /£.»)/#  -  (/#  + 


Next  considering  (3.8),  we  have 

<  KNfN,  CNfN  >RN=  N  {2fy  -  f»)f»  +  2  N2(f%  -  /^_j)2 


Ar- 1 


+*  E  *  ((A"  -  AM  -  (At.  -  A"))  ((A"  -  A"-,)  +  (/, 


Jt+1 


*  =  2 
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In  a  manner  similar  to  that  above,  the  sum  £*-2  can  be  rewritten  as 
N-l 


N  £  [(*  - 1  )(A"  -  A"-,)’  -  *(/,"♦.  -  A")3  +  (A*  -  A"-,)3] 

k=2 

=  -JV(JV  -  i)(/£  -  /£_,)’  +  w(AK  -  /")*  +  Jv  £(a"  -  A"-,)3- 

k=2 

Using  the  above  equality,  we  have 

<  K’ti\c”r  >**=  N  ((A")3  +  £(A*  -  A"-.)3)  +  *3(A?  -  //?-.)’• 

Using  Lemma  3.1,  we  obtain  (3.8). 

By  the  definition  of  the  matrix  BN  we  have  for  any  /N,g'v  6  RN , 

<  BNfN,fN  >rN=  a{f%)\  <  BNfN,CNgN  >R»  =  2 aNf»(g%  -  g^). 

Combining  the  results  in  Lemma  3.3,  we  can  obtain  Condition  2.3. 

Lemma  3.4  Condition  S.S  holds. 

Proof:  Consider 


□ 


<MY,/N>^  +  <MY,CYv  =  ~~{fs  +  /n-i)2  +  /jv_i) 


<  4(/")2  =  -  <  BNfN,fN  >Ry  . 

a 


Also,  we  have 


<  /fVY  >**  -  <  K"g\CNgN  >R»  -  <  BN fN >CN gN  >RN 

—  ~^7(9n  ~~  9n-i)7  ~  2&N Jn{9n  ~  Sn-i) 

<a'(f»y  =  a<BsfN,fN  >R». 

□ 

We  can  thus  invoke  Theorem  2.1  to  conclude  that  the  method  described 
above  preserves  uniformly  exponential  stability  of  the  solutions  of  equations  (3.1)  - 
(3.3). 
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In  the  case  of  a  two-dimensional  wave  equation  on  a  rectangular  domain 
ft  =  {(*>y)  :  0  <x<l,0<y<l},we  may  consider  the  following  special  case  of 
(1.1)-(1.3): 

d 2  d 2  d 2 

(3.9)  3^u(f,x,y)  =  ~u(t,x,y)  -)-  —  u(t,x,y),  (x,y)  £  Cl,t  >  0, 


(3.9) 

(3.10) 


(3.10)  u(t,x,y)  =  0,  (x,  y)  6  T0,  t  >  0, 

Q 

(3.11)  dr,u(t,x,y )  =  -a—  u(t,x,y),  (x,y)6T1,t>  0, 

where  F0  =  {(x,y)  €  ft,x  =  0,  or  y  =  0}  and  Tj  =  n/T0.  The  mixed  finite  element 
method  is  defined  as  follows.  For  each  given  integer  N,  let  (x,-,y3)  be  the  grid  points 
defined  by  x,  =  i/7V,  y3  —  j/N,  for  i,j  =  0,1,...,  N.  For  each  pair  of  indices  (i,j) 
with  1  <  t,  j  <  N,  a  neighborhood  of  the  grid  point  (x,-,y3)  is  defined  as 

D"  =  {(^,y)  e  ft  :  |x  -  X,|  <jj,\y-  y3|  <  |x  -  y  -  xt-  +  y,-|  < 

Two  families  of  functions  are  defined  by 

j.Nt„  _  /  1  -  iVmax{|x  -  x,|,|y  ~  y3|,|x- y  -  x,  -  y3|},  (i,y)eDj, 


<M*>y)  = 


v  ’  '  (0,  otherwise, 

ibN(x  v)  =  |  lij'  (x>y)eAy» 

,J  {  0,  otherwise, 

where  the  constants  7 tJ  are  chosen  such  that  we  have  the  equality 

f^d"(x,y)dxdy  =  j^ip£(x,y)dxdy. 

We  then  choose  =  span{^}£y_.,,  =  span{^}(vJ=r  We  prefer, 

for  simplicity  of  presentation,  to  denote  an  element  of  an  N 2  x  N2  matrix  AN  by 
Using  this  notation,  the  matrices  ,KN ,  and  BN  axe  defined  by 


Km  = 

Ju  'JJ"{x,y)'Jj%{x,y)dxdy, 

TS  N  _ 

K  (•>).(«)  - 

fn  V<t>”{x..y)  ■  V<f>%(x,y)dxdy, 

dN  _ 

[  a<f>?J(x,y)4>%{x,y)dT1. 

J  Pi 

Although  our  numerical  results  indicate  that  this  method  preserves  uniformly  expo¬ 
nential  stability  of  the  solutions  of  (3.9)-(3.1l),  we  are,  to  date,  unable  to  confirm  it 
analytically.  We  have  attempted  to  use  arguments  analogous  to  the  one-dimensional 
case,  but  the  number  of  terms  to  compute  and  the  associated  tedium  quickly  become 
overwhelming. 
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4  Polynomial  based  Galerkin  approximation 


In  this  section,  we  present  a  stability  analysis  for  a  general  class  of  Galerkin  type 
approximations  of  the  wave  equation  in  Rn.  Several  technical  assumptions  on  the 
domain  f!  and  the  approximation  spaces  H? ,  H £  are  needed  for  exponential  stabil¬ 
ity  of  the  approximate  solutions.  In  particular,  hypercubic  domains  in  Rn  with  all 
reflecting  sides  sharing  at  least  one  common  corner  satisfy  these  conditions.  The  re¬ 
quirements  on  the  approximation  spaces  can  be  satisfied  if  we  take  Hfi  =  H =  HN , 
and  HN  is  a  carefully  defined  subspace  spanned  by  the  tensor  products  of  polyno¬ 
mial  functions. 

Since  the  arguments  developed  here  can  also  be  useful  for  the  analysis  of 
other  approximation  methods,  the  assumptions  are  made  in  a  general  form.  We 
will  indicate  in  each  case  how  polynomial  based  Galerkin  methods  satisfy  these 
assumptions. 

Now  consider  the  wave  equation  given  by  (1.1)-(1.3)  where  we  make  the 
following  assumptions  on  the  domain  Cl. 


Assumption  4.1  There  exists  a  constant  r  >  0  such  that  for  any  e  >  0,  there 
exists  a  vector  field  £(■)  6  C4{Cl\Rn),  with  the  properties: 

(Vl)  The  matrix  L{x)  defined  by 


+  «,,(*)) 


satisfies  L(x)  -  /  >  0  on  the  domain  Cl; 

(V2)  For  x  £  dCl,  the  outward  normal  unit  vector  rj(x)  satisfies 


£(x)  ■  rj(x)  <  0,  forxeTo,  £(x)  ■  r}[x)  >  r,  for  x  G.  T \\ 
(VS)  The  following  inequalities  hold 

I  it  ka\  —  g  I  T,  km  I  <  f- 

«.i= i  =  l 


where 


d 2  a3 

a-  M1)'  k{?i  ~ 


dxxdxj 


dx,d2Xj 
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We  note  that  (Vl)-(V3)  are  exactly  the  same  conditions  as  these  required 
in  [Cj  where  the  author  treats  the  exponential  decay  of  the  solution  of  the  wave 
equation.  For  a  given  domain  fl,  if  there  exists  a  point  x0  £  Rn  such  that  (x  - 
x0)  •  r?(x)  <  0  on  To  and  (x  —  x0)  •  r?(x)  >  r,  on  IT,  by  taking  £(x)  =  x  -  Xo,  the 
conditions  (Vl)-(V3)  are  clearly  satisfied.  In  particular,  consider  the  hypercube  f 1 
given  by  fl  =  {x  €  Rn,0  <  x0  <  1};  if  T0  is  the  union  of  a  collection  of  sides  that 
share  at  least  one  common  point  x0,  then  one  can  verify  that  £(x)  =  x  —  xo  satisfies 
(Vl)-(V3).  In  fact,  in  the  case  of  a  hypercube,  if  r0  is  the  union  of  a  collection  of 
sides  without  any  common  point,  it  can  be  shown  that  the  solutions  of  the  partial 
differential  equation  (l.l)-(l.3)  are  not  exponentially  stable. 

Now  consider  subspaces  H k  =  HN  C  H2(Cl)  Pi  H^o(Cl)  given  by 

II N  =  spanl^^J^-j,  where  the  <fk  are  a  general  family  of  approximation  elements. 
The  matrices  MN  ,K*  ,  BjV  of  Section  2  are  defined  by 

=  ftf  (x)4>H*)dz, 

K$  =  fv^lx)-V^(z)dx 


Before  we  introduce  the  multiplier  matrix  CN ,  we  present  a  few  useful  special 
properties  of  the  matrices  MN ,  KN ,  BN . 

Lemma  4.1  Consider  a  sequence  { ( ' ) } jt*L i  °f  functions  in  HN  and  let  V*  be  the 
matrix  defined  by 

Kf  =  ft?  (x)6?  (z)dx. 

Then,  the  matrix  DN  =  (V^ )T  [MN)~i  BN  can  also  be  given  by 

D*  =  L  ae>  (x)4>? {x)dx. 

JY  i 


Proof:  Since  the  6 k  (-)’s  are  elements  of  /fA,  there  exist  constants  6fk  such  that 

j= i 

Let  QN  be  the  matrix  with  entries  9^};  then  we  have  =  MNQN.  In  fact. 


k?  =  e<  f  =  f; «. 

1=1  1=1 
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Therefore,  we  obtain  DN  =  ( QN)TBN  and  hence 


Dn  = 


/  E<(*)<(*)^=  f  O0?{z)tf{x)dx. 

1=  1  'r‘i=l  ■'*» 


□ 


Lemma  4.2  Let  be  the  same  matrix  as  defined  above.  Then,  the  matrix  FN  = 
(VN)T (MN)~lKN  is  also  given  by 


F; 


jv 


-^<j>?(x)dx  -  ^  0*(x)A^f  (x)dx. 


Proof:  We  only  need  to  consider  =  [QN)TKN  where  the  matrix  0W  is  the  same 
as  before.  We  have 


f;]  =  f  E^V0r-v^(x)rf* 

yn  i=i 

=  [  (x)dx  -  [  (x)da 

■'Ti  C77J  Jn  j_ J 

Therefore,  we  obtain 

FS  =  ^  {x)dx  -  0?{x)  ■  A<f>? ( x)dx . 


□ 

Now,  by  imitating  the  multiplier  used  in  Section  2  for  Q(t),  we  define  the 
matrix  WN  —  (CN)T MN  as 

k"  =  fm"  ■  +  <£x,  - 

■/n  k=i 

where  £N^-)  is  a  vector  field  which  satisfies  Assumption  4.1.  In  order  to  be  able  to 
use  the  previous  matrix  equalities,  we  need  further  technical  assumptions  on  £Ar. 


Assumption  4.2  There  exist  constants  r  >  0  and  M  >  0  independent  of  N  such 
that  there  exists  a  vector  field  £N (•)  £  C*(Cl;Rn)  which  satisfies  Assumption  4-1  and 
furthermore 

2 lN  •  V4Z  +  XX, <  G  HN,  with  ||£"(x)|U„  <  AT,  !XX(*)I  <  M, 

i=i  i=i 

for  all  x  €  0  and  /c  =  1 ,  •  •  • ,  TV . 
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The  above  assumption  is  satisfied  if  there  exists  x0  £  Rn  such  that  ts  (x)  = 
x  —  xo  satisfies  Assumption  4.1  and  if  HN  is  defined  by 

HN  =  |  /(*!,..., *„)  =  f[Pk{xk)  6  #r0(fi),  where  1 

(.  pk(-)  are  polynomials  of  degree  less  than  or  equal  to  mi.  J 

In  the  case  of  a  hypercube,  one  can  easily  construct  HN  by  taking  p*  £  Pmt(0,  l), 
where  Pmjt(0,l)  are  polynomials  of  degree  less  than  equal  to  m*  which  satisfy  the 
appropriate  Dirichlet  boundary  condition. 

The  following  lemma  will  establish  Condition  2.2. 

Lemma  4.3  Under  Assumption  4-1,  there  exists  constants  Cj,C2  depending  only  on 
the  domain  H  such  that 

I<mY,cVvI  =  |  <  wW*  >*»  I 

<  c1<KNgN,gN  >RN  +c2<MNfN,fN  >Rs, 

for  all  fN,gN  £  RN . 

Proof:  Let  us  associate  each  vector  fN  £  RN  with  an  element  fN  (•)  £  HN  given  by 
fN(x)  =  I2k=i  fk  4>k  (x)-  ^  is  easy  to  verify  that  for  any  vectors  fN ,gN  £  RN 

<KNg",  =  /nl|V9"(I)||i.<ix  =  ||9''(.)||;, 

Similarly,  we  have 

<  MNfN,CNgN  >RN  =  <  WNfs,gN  >Rs 
=  f  2/"(x)  (*"(*)  •  Vg«(x))  +  (£  (*)  -  1  j  /"(*)<,"(*)  iz. 

From  standard  results,  there  exists  a  constant  c  depending  only  on  Q  such  that  for 
all  u  £  Z/po(n),  ||u|||2(nj  <  cjjtzjjf .  By  Assumption  4.2,  we  obtain 

I<a/'7",cV'  >»,|  <  5ii/"()|ii=(„) 

+  (j  £<&(•) -1  +»il<K(-)IU-(n)|ll»''(-)ll? 

v  *=i  L°°(n)  J 

=  c2  <  MNfN,fN  >RN  +d  <  KNgN,gN  >RN, 

where  c2  =  3/2, C!  =  ™||£Ar(‘)iUo°(n)  +  c/2||  53"=  i  ^(')  _  l||L°°(n)-  □ 

The  following  two  lemmas  are  useful  in  establishing  Condition  2.3. 
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Lemma  4.4  Let  Assumption  4-2  hold.  For  any  fN  £  RN ,  let  fN  (•)  fee  the  function 
in  Hn  given  by  fN{x)  =  E?=1  Then 

<  MNfN,CsfN  >rN  =  <WNfN,fN  >rN 

=  -  <  MN fN  JN  >RS  +  f  fN(x)2CN(x)  ■  T](x )dx. 

■T, 


Proof:  Let  GN  =  Ws  +  M*;  then 

Gi?  =  fn  {2(<"  •  +  (E  <t)*f  (x)*/(x) }  dx. 

Then,  we  have 


<GNfNJN>RN 

N 


=  fa  E  {2(^  ■  W?//X(*)/f  +  (EmW)  (*)/* |  dx 

=  fn  {2^"  ■  V/N(x))/*(x)  +  E<tW/"(x)2J  dx 
=  J^d\v(fN(x)2tN  (x))dx 
=  Jr  fN(x)HN{x)-r]{x)dx, 


because  fN{x)  ~  0  on  r0. 


□ 


Let  us  define  the  matrix  L^x)  by  (recall  condition  (Vl)  of  Assumption  4.1) 


dij 


3x,  ) 


Lemma  4.5  Under  Assumption  4-2,  for  any  gN  €  /ef  g^(-)  fee  f/ie  function  in 
HN  given  by  gN  (x)  =  E*Li  9k<f>kix)-  Then  we  have 

<  KNgN,CNgN  >RN  -  <  jfV,^  >** 

=  2  /n(Vs")r(£."(x)  -  /)V3"(x)dx  -  i 

+  f  IVT*"  ■  .jdx. 

■'r‘  x  ,;=1  •'r, 
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Proof:  We  have  <  >*»  =  <  >„.  By  ,aking 


0;VM  =  2«"  •  v<*"W  +  ^  <£*(*)  - 1  ]  w, 


we  have 


where  K*  is  as  defined  in  Lemmas  4.1  and  4.2.  By  Lemma  4.2,  we  have 

where  the  matrix  is  given  by 

F*  =  Lx  ^dx  ~  fu  **(*)  ■  M?(x)dx. 

As  a  consequence,  we  obtain 

(4.i)  <*V,cV  >„„ 

=  -/„(«"  ■  V5»  +  (£«"  _  1)9»(I))A,»(I)<(I 

*=  1 

+/r,  (2<" '  *»"(*)  +  (E*  -  A, 

Consider  first  the  term 

/"  =  /n  («»  .  Vp"(x)  +  (g<»  -  !),*(*)) 

For  all  functions  4>  6  H2(U) 

2Af,((»  ■  TV)  +  (gf*  j  ^ 


div  |2V*(<"  ■  VO)  -  IVOIV  +  e^)V4, 

1=1 

-2(vo)ri"vo  -  22  €  *—■ 


Therefore,  we  obtain 


IN  =  -fji^gN(z)-LNVgN(x)+f2l^gNj^+gNAgN(x)^ 


dx 


+ 


k=  1 


2  (^-V,*) +  (£<,)* 

-  [  |VS"|V  •  D<far, 

Jra 


N  \N  \  in„N\2cN 


dv 


-\VgN\2r  -r, 


dx 


where  we  have  used  Assumption  4.2  to  observe  that  2{tN  •  V gN)  +  ££=1  tk,k9N  ip  in 
Hn  and  hence  vanishes  on  To-  Since  gN  vanishes  on  To,  we  have  that  also  2 [LN -VgN) 
must  vanish  on  To-  Moreover,  since  gN(x)  =  0  on  To,  V gN (x)  =  ±|V(?n|t7  on  Tc 
and  hence  |VgjV|2Tv  •  rj  =  1 1 Vg,v \tN  ■  VgN  =  0  on  To-  As  a  consequence,  we  obtain 


IN  = 


2  <  LNVg\VgN  >rN  +  Y:  ^°§y/  +  *****  f  dX 


«j=i 


L 


2(£^-v^)  +  (E*)^ 


dg 


N 


'VgN\HN-rt 


Jk=l 


dr\ 


dx. 


Integrating  the  terms 


N 


0  =  1  OX3 


N\2 


by  parts  once,  we  obtain 

/"  =  -  /  { 2  <  LNVgN ,  V9"  -|V/T  -  (  £ 

jn  0  =  1  2 

2£‘V  ■  V/’  +  (Y  £&  -  1  )g> 

k=  1 


I  drj 

(gNy 


dx 


dx. 


0=1 


Finally  combining  the  above  result  with  (4.1),  we  obtain 
<  KNgN,CNgN  >rN 

=  f  \ 2  < 

;n  ,^i  2 

+  /rj||VS,''p£"  + 


<hr 
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This  least  equality  is  equivalent  to 


<  KNgN,CNgN  >RN  -  <  KNgs,gN  >R» 
f  2(<  LNVgN  ,VgN  >„„  - 1  VS~|!)  dr 

-l  Jjt  +  L  |VT<"  '  nd* 

l  jq  ij=1  j r\ 


By  using  Lemma  4.1  and  the  same  notation  as  for  Lemma  4.4,  4.5,  we  find 

<b"/n,  =  /ri«U''.vs'  +  £(;i-i)r'j/'i, 

<  BNfN,fN  >rn  =  [  a(!»fdz. 

Jr, 

Lemma  4.6  Under  Assumption  4-2,  Condition  2.8  holds. 

Proof:  (i)  By  the  boundedness  of  iN  and  Lemma  4.4,  (Si)  holds,  (ii)  By  taking  e 
sufficiently  small  in  Assumption  4.1  and  4.2,  we  have 

fTi  jlv/r(*)|2^Ar(*)  ^{x)  +  Ci^)^'V(2)2|  dx  >  0, 


and 


\  fit  <  « l  ir7g>1(x)\,dx 


for  5  arbitarily  small.  Therefore  (S3)  following  Condition  2.3  holds,  (iii)  By  the 
Schwartz  inequality,  <  BN fN ,CN gN  >rn  can  be  bounded  by 

|<5Y,cVvl<-/  \f"(x)\2dx  +  c2o  [  \VgN(x)\2dT 

V  J r,  Jr, 

where  the  constants  cuc2,  and  v  are  independent  of  N  and  u  can  be  arbitrarily 
chosen.  Therefore,  (S2)  of  Condition  2.3  holds. 

As  a  consequence  of  Lemma  4.3  and  4.6,  on  the  domain  fl  and  the  sub¬ 
space  HN ,  we  conclude  that  under  Assumption  4.2  these  approximation  methods 
uniformly  preserve  the  exponential  stability  of  the  weakly  damped  wave  equation. 


5  Numerical  studies  of  approximation  methods 


As  has  been  demonstrated  in  the  previous  three  sections,  establishing  a  uniform 
decay  rate  estimate  for  a  given  approximation  method  can  be  complex  and  tedious. 
In  fact  the  general  approach  outlined  in  Section  2  offers  only  a  broad  direction 
toward  the  selection  of  the  appropriate  multiplier  QN .  On  the  other  hand,  for  any 
given  approximation  method,  by  computing  the  eigenvalues  of  the  matrices  AN  in 
(1.5),  one  can  observe  the  trends  in  the  location  of  the  eigenvalues  as  N  increasess. 
In  several  cases  presented  below,  the  numerical  results  clearly  indicate  that  some  of 
the  eigenvalues  of  AN  are  approaching  the  imaginary  axis  as  N  increases,  and  it  is 
therefore  unlikely  that  a  uniform  decay  rate  can  be  preserved. 

The  example  used  for  demonstration  in  this  section  is  the  two-dimensional 
wave  equation  (3.9) - (3. 1 1 )  where  the  parameter  a  is  taken  to  be  1.  The  construction 
of  the  matrix  .4‘N  is  presented  for  each  approximation  method.  The  eigenvalues  are 
then  computed  using  the  subroutine  F02AFF  of  the  FORTRAN  software  library 
NAG.  The  results  were  compared  with  the  results  using  the  IMSL  library,  and  no 
visible  difference  can  be  observed.  All  computations  were  carried  out  on  a  IBM 
3081  computer. 

5.1  Polynomial  based  Galerkin  approximation  methods 

Let  {Q/fc(-)}£to  be  the  Legendre  polynomials  of  degree  less  than  or  equal  to  N  on 
the  interval  '-1,1'.  We  define  a  sequence  of  polynomials  P*  by 

=  <?*(*) +  (-l)*Qo(s) 


for  k  = 

It  is  easy  to  see  that  P,N  —  span{P*}£Li  is  the  set  of  polynomials  of  degree 
less  than  or  equal  to  .V  vanishing  at  s  —  —1.  Let  us  define  matrices  A/.  ,  as 
follows 

I  AC’.,  =  f\' P,(s)PAs)ds, 

The  matrices  M .v ,  K*  can  be  easily  computed  using  the  recursive  properties  and 
the  orthogonality  of  the  Legendre  polynomials.  Now  let  dP(x,y)  =  dP|(x,y)  = 
P,(x)P;(y).  We  assume  a  scaling  of  the  domain  Q  onto  [  —  1,1]  x  [—1,1]  is  made. 


The  matrices  Ms ,KN ,  and  BN  defined  in  the  previous  section  can  be  computed  as 
follows: 


iA"'|W).(*0  = 

[B"]W>.(«)  =  a|M,A']jl/>,(l)C,(l)  +  a/>,(l)n(l)(M,"!„. 

The  matrix  AN  is  given  by 

As_(  0  /  ^ 

\  ~(Mn)~1Kn  ~(Mn)-'Bn  )  ■ 

The  dimension  of  AN  is  N 2.  For  N  =  3,4,6,  ,9,  the  locations  of  eigenvalues  are 
displayed  in  Figure  5.1.  We  note  that  as  would  be  predicted  from  the  analysis  of 
Section  4,  a  uniform  margin  between  the  eigenvalues  of  AN  and  the  imaginary  axis 
is  maintained  for  all  N.  In  the  Table  5.1,  we  list  the  margin  of  stability  for  each  N. 


Table  5.1:  Margin  between  the  eigenvalues  of  the  matrix  AN  and  the  imaginary 
axis. 


N 

max{ReA,A  6  a(AN)} 

4 

-0.6230 

5 

-0.6215 

6 

-0.6232 

7 

-0.6097 

8 

-0.5815 

9 

-0.5628 

5.2  Polynomial  spline  based  Galerkin  approximation  methods 

For  N  a  given  integer,  the  interval  [0,  l]  is  divided  into  equal  sized  subintervals 
\xk-i,Xk\  with  xk  —  k/N  for  A:  =  1  Let  BN,m  be  the  set  of  polynomial 

splines  of  order  m  corresponding  to  the  grid  {x*}  that  vanish  at  x  =  0.  Consider 
a  basis  for  BN'm  on  ft  =  [0,  l]  x  [0,  lj.  The  matrices  M*,K?  are 

defined  as 

[<]„  =  f1  B*'m(x)B”‘m(z)dz,  i,j  =  l,...,N  +  m-  1, 

Jo 
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The  matrices  M*,K?  can  be  easily  evaluated.  We  define  4>*(z,y)  =  t P-j(x,y)  = 
B?'m(x)B^'m(y).  The  the  matrices  MN,KN  and  BN  are  defined  as 

[ma’](.j),(*/)  =  <f>*  {x,  y)4>%(x,  y)dxdy, 

\KN)[tjUkl]  =  f  V<f>^{x,y)  ■  V<f>%{x,y)dxdy 
J  a 

[BN}[tj),(ki)  =  ^a<t>”j[x,y)<t>Z{x,y)dr1. 

Again,  the  matrices  ,KN ,  and  BN  can  be  computed  using  and  K* .  Since 
only  the  first  derivative  of  the  spline  function  is  required,  either  linear  spline  or 
cubic  spline  functions  can  be  used.  In  Figures  5.2  and  5.3,  the  locations  of  the 
eigenvalues  of  the  matrix  AN  using  linear  and  cubic  spline  functions,  respectively, 
are  depicted.  We  note  the  similarity  between  the  locations  of  the  eigenvalues  with 
small  imaginary  part  and  the  locations  of  the  eigenvalues  of  the  matrix  AN  using 
polynomials  depicted  in  Figure  5.1.  However,  the  eigenvalues  with  large  imaginary 
part  tend  toward  the  imaginary  axis.  Another  interesting  observation  is  that  un¬ 
like  the  eigenvalues  of  using  polynomials  and  cubic  splines,  the  eigenvalues  of 
AN  using  linear  splines  do  not  have  exceedingly  large  negative  real  parts;  further 
discussion  concerning  this  observation  will  be  given  in  the  final  section. 


5.3  Finite  element  method 

The  classical  finite  element  method  requires  that  we  subdivide  the  domain  Cl  into 
triangles  where  functions  linear  on  each  triangle  are  used  for  the  approximation. 
More  precisely,  for  any  1  <  i  <  N  and  1  <J<  N,  a  function  <f>^  is  defined  as 

wv,  _  /  l-Armax{|z-x,|,|y-yJ|,|x-y-x,+yi|},  (z,y)el$, 

'jo,  otherwise, 

where  the  support  D -j  is  defined  by 

Dii  =  {(z.y)  e  n:\x-  x,|  <  -^,|y  -  y>|  <  y-  x,  +  y,-|  <  -^} 

with  x,  =  i/N,  y3  =  j/N.  We  choose  the  functions  rp^(x,y)  =  <f>ij{x,y)  and  the 
matrices  MN,KN,  and  BN  are  defined  as  in  the  other  Galerkin  methods. 
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The  computation  of  the  matrices  MN,KN,  and  BN  uses  finite  element  tech¬ 
niques  based  on  elementary  matrices. 

The  location  of  the  eigenvalues  are  depicted  in  Figure  5.4.  We  note  that 
some  eigenvalues  with  large  imaginary  part  tend  toward  the  imaginary  axis  which 
suggest  that  no  uniform  decay  rate  of  solutions  is  preserved.  We  also  observe  that 
only  a  few  eigenvalues  have  large  negative  real  part. 

5.4  Mixed  finite  element  method. 

The  mixed  finite  element  method  presented  in  Section  3  was  also  implemented.  In 
Figure  5.5,  the  location  of  eigenvalues  are  depicted. 


Figure  5.5:  Location  of  the  eigenvalues  of  the  matrix  AN  for  the  mixed  finite  element 
method. 


Although  a  uniform  margin  between  the  eigenvalues  of  the  matrix  AN  and 
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the  imaginary  axis  is  maintained,  the  locations  of  the  eigenvalues  are  dramatically 
different  from  the  previous  approximation  schemes  at  lower  values  of  N .  This  might 
suggest  a  slower  convergence  rate  for  the  approximation  scheme. 

In  Table  5.2,  the  size  of  the  margin  between  the  eigenvalues  of  AN  and  the 
imaginary  axis  is  reported. 


Table  5.2:  Margin  between  the  eigenvalues  of  the  matrix  AN  and  the  imaginary  axis 
for  the  mixed  finite-element  method. 


N 

max{ReA,A  E  o(AN )} 

4 

-0.4689 

5 

-0.4596 

6 

-0.4552 

7 

-0.4527 

8 

-0.4512 

9 

-0.4969 

10 

-0.4496 

15 

-0.4632 

5.5  Finite  difference  approximation  method 

Unlike  the  other  approximation  methods  presented  so  far,  in  the  finite  difference 
method  one  attempts  to  approximate  the  solution  of  the  wave  equation  at  the  grid 
points  directly.  In  fact,  for  any  given  integer  N,  let  u* (t) ,  (t)  be  the  approxi¬ 

mation  of  the  solution  values  u(t,  z,,yj),v(t,  z,y:)  =  ut(t,  x,,  y;)  of  (3.9)-(3.1l).  The 
wave  equation  is  approximated  by 

(5.D  jt<m 

(5.2)  !<(«) 

for  i  =  2, . . . ,  N  —  1  and  j  ~  2, . . . ,  N  —  1. 

At  the  boundary  z  =  0  or  y  =  0,  by  taking  u$:  =  u^0  =  0,  the  equation  (5.2) 
still  holds.  However,  at  the  boundary  i  =  1  or  y  =  1,  the  second  derivatives  are 


=  t&(0. 


n2  (~4uJ>(0  +  u£ij(0  +  «?+.j(0  +  “S-i  +  «5+i(0)  i 


3-1 


approximated  by  one  sided  finite  differences.  In  particular,  we  take 

=  JV>  (-2<„(i)  +  «£,,„(<)  +  ««.,„(!)) 

+  N  (av?N(t)  -  N(u?N(t)  — 

=  N 1  (-*<,(«)  +  <,_,(()  +  «S„>.(0) 

+N  (aujlrj(0  _  ^(ujv,j(0  ~  ujv-i,j(0)) 

for  t  =  l,...,jV  -  l,j  =  1,...,JV  -  1.  At  the  corner  (1,1),  we  take 

=  N  (at;Ar,.v(0  ~  N(ugN(t)  ~  u$,N-l(0)) 

+7V  (a<N(t)  -  N(u^(t)  -  ug_hN(t)))  . 

The  locations  of  the  eigenvalues  of  the  matrix  AN  are  depicted  in  Figure  5.6. 

It  was  a  surprise  to  us  that  the  eigenvalues  of  AN  for  this  finite  difference 
approximation  are  so  different  from  those  of  the  other  schemes.  Many  variations  of 
the  finite  difference  methods  exist;  we  believe  that  an  adjustment  of  the  method  used 
here  may  well  produce  different  locations  for  eigenvalues.  However,  our  experience 
with  the  above  standard  finite  difference  approximation  method  suggests  that  no 
uniform  margin  between  the  eigenvalues  of  AN  and  imaginary  axis  is  maintained 
for  all  N. 
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6  Concluding  remarks 


Even  in  view  of  the  analysis  for  the  methods  given  in  Sections  3  and  4,  one  might 
ask  why  these  two  methods  succeed  while  several  other  popular  approximation 
schemes  fail.  One  possible  explanation  is  the  following.  As  approximation  methods 
of  a  second  order  system,  many  approximation  methods  approximate  the  second 
component  v  which  has  the  same  smoothness  as  Vu  theoretically  by  functions  with 
equal  order  of  smoothness  as  the  approximation  for  u.  In  another  words,  if  E 
and  Vq  E  H £  are  taken  to  be  the  initial  values  of  the  equations  (l.l)-(1.3),  the 
solution  t>  is  much  less  smooth  than  elements  in  .  This  is  not  the  case  in  the 
mixed  finite  element  method  and  the  polynomial  based  Galerkin  methods. 

We  should  also  comment  on  the  implications  of  the  present  results  with  re¬ 
spect  to  the  boundary  control  problem.  If  the  boundary  condition  (1.3)  is  considered 
as  a  boundary  feedback  control,  with  an  appropriate  approximation  of  the  input 
operator,  the  results  reported  here  can  be  used  to  construct  finite  dimensional  ap¬ 
proximate  control  systems  that  are  uniformly  exponentially  stabilizable.  However, 
this  does  not  guarantee  that  the  infinite  dimensional  system  under  the  approxi¬ 
mate  control  is  exponentially  stable.  The  inherent  sensitivity  of  the  conservative 
wave  equation  with  respect  to  perturbations  of  the  boundary  conditions  remains 
one  major  drawback  of  this  type  of  model. 

Finally,  we  observe  that  in  both  polynomial  based  Galerkin  methods  and 
the  mixed  finite  element  method,  a  large  number  of  extraneous  eigenvalues  with 
large  negative  real  parts  are  introduced.  Therefore,  the  resulting  finite  dimensional 
ordinary  differential  equation  generated  by  these  methods  can  be  very  stiff.  This  can 
cause  problems  in  integration  algorithms.  Since  our  main  concern  is  to  solve  linear 
quadratic  control  problems  using  these  approximation  schemes,  this  observation 
should  not  substantially  diminish  the  attractiveness  of  these  two  methods. 
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